class BeringWaterCol extends Slice {
  
  float[] PARforcing, forcingTime, forcingZ, forcingKfrac, forcingH, iceForcing;
  int[] forcingKbefore;
  float[][] Tforcing, kappaForcing;
  
  BeringWaterCol(int K, float Htot, float dz0, String paramFilename, String copepodName, float tStart, float tEnd, ForcingSet forcing, int forcingIndex) {
    name = "Bering Sea water column";
    shortname = "LowLaMB";
    
    float ddz = 2 * (Htot - K * dz0) / K / (K-1);
    float[] dz = new float[K];
    for (int k=0; k<K; k++) dz[k] = dz0 + (K-1-k)*ddz;
    createDomain(new float[] {1}, dz);
    
    linkToEcosystem(new BeringEcosystem(paramFilename, copepodName));
    ((BeringEcosystem)eco).env = this;

    Statevar T = addVar("temperature","temp","°C",Ncells);
    T.advects = false;
    T.mixes = false;
    
    allocateStorageForPhysics();
    
    addParam("surface PAR","PAR0","W/m2",150,0,300).hidden = true;
    addParam("fractional ice cover", "ice", "", 0, 0, 1);
    addSinking("SD", getParamVal("wsink_SD"), 1);
    addSinking("LD", getParamVal("wsink_LD"), 1);
    
    // assemble forcing
    forcingTime = forcing.time();
    forcingZ = forcing.z();
    forcingH = forcing.one_H_case(forcingIndex);
    for (int i=0; i<forcingH.length; i++) forcingH[i] = abs(forcingH[i]);
    PARforcing = forcing.one_PAR_case(forcingIndex);
    Tforcing = forcing.one_T_case(forcingIndex);
    kappaForcing = forcing.one_kappa_case(forcingIndex);
    iceForcing = forcing.one_ice_case(forcingIndex);
    // save info to make it easy to interpolate in depth later
    forcingKbefore = new int[forcingZ.length];
    forcingKfrac = new float[forcingZ.length];
    for (int k=0; k<K; k++) {
      float zk = constrain(z_rho[k], min(forcingZ), max(forcingZ));
      int[] fk2 = findBeforeAfter(forcingZ, zk);
      forcingKbefore[k] = fk2[0];
      forcingKfrac[k] = constrain((zk-forcingZ[fk2[0]]) / (forcingZ[fk2[1]]-forcingZ[fk2[0]]), 0, 1); // fractional index 
    }
    
    // set i.c.s
    float Hcurrent = interp1(forcingTime, forcingH, tStart);
    for (int k=0; k<K; k++) {
      getVar("NO3").initial[k] = ((BeringEcosystem)eco).NO3bc(z_rho[k], Hcurrent);
    }
    getVar("NH4").setInitial(0.0);
    getVar("SP").setInitial(getParam("ic_SP").initial);
    getVar("LP").setInitial(getParam("ic_LP").initial);
    getVar("MZ").setInitial(getParam("ic_MZ").initial);
    String[] stages = ((BeringEcosystem)eco).stages;
    for (int i=0; i<stages.length; i++) getVar(stages[i]).setInitial(getParam("ic_"+stages[i]).initial);
    getVar("Eu").setInitial(getParam("ic_Eu").initial);
    getVar("SD").setInitial(0.0);
    getVar("LD").setInitial(0.0);
    
    time.setInitial(tStart);
    stopTime = tEnd;
    startArchive(0.5);
    for (int i=0; i<fluxes.length; i++) {
      String theName = fluxes[i].shortname;
      if (theName.length() > 5 && theName.substring(0,5).equals("hDiff")) archive.exclude(fluxes[i].shortname);
      if (theName.length() > 7 && theName.substring(0,7).equals("lateral")) archive.exclude(fluxes[i].shortname);
      if (theName.length() > 4 && theName.substring(0,4).equals("hAdv")) archive.exclude(fluxes[i].shortname);
      if (theName.length() > 4 && theName.substring(0,4).equals("vAdv")) archive.exclude(fluxes[i].shortname);
    }
    // save the raw forcing to the archive
    archive.addNameValuePair("forcingTime",forcingTime);
    archive.addNameValuePair("forcingZ",forcingZ);
    archive.addNameValuePair("forcingH",forcingH);
    archive.addNameValuePair("Tforcing",Tforcing);
    archive.addNameValuePair("kappaForcing",kappaForcing);
  }


  float icecover() {
    return interp1(forcingTime, iceForcing, time.current[0]);
  }
  float ddt_icecover() {
    return interp1(forcingTime, iceForcing, time.current[0]+0.5) - interp1(forcingTime, iceForcing, time.current[0]-0.5);
  }
  
  void interpForcingInTime() {
    getParam("PAR0").current = interp1(forcingTime, PARforcing, time.current[0]);
    getParam("ice").current = icecover();
    float minKv = pow(10,-5.5);
    float maxKv = pow(10,-2); // for FTCS, stability limit is 1/2 * grid spacing^2 / timestep
    // adjust water column depth
    float Hcurrent = interp1(forcingTime, forcingH, time.current[0]);
    int kbot_old = kbot[0];
    for (int k=0; k<K; k++) if (z_rho[k] < -Hcurrent) kbot[0] = k;
    kbotu[0] = kbot[0];
    kbotu[1] = kbot[0];
    // actually interpolate forcing
    for (int m=0; m<M; m++) {
      for (int k=0; k<K; k++) {
        int ind = km2ind(k,m);        
        float[] yy = new float[2];
        yy[0] = interp1(forcingTime, Tforcing[forcingKbefore[k]], time.current[0]);
        yy[1] = interp1(forcingTime, Tforcing[forcingKbefore[k]+1], time.current[0]);
        getVar("temp").current[ind] = (1-forcingKfrac[k]) * yy[0] + forcingKfrac[k] * yy[1];
        yy[0] = interp1(forcingTime, kappaForcing[forcingKbefore[k]], time.current[0]);
        yy[1] = interp1(forcingTime, kappaForcing[forcingKbefore[k]+1], time.current[0]);
        float kvk = (1-forcingKfrac[k]) * yy[0] + forcingKfrac[k] * yy[1];
        kappa[k][m] = constrain(kvk, minKv, maxKv) * 86400;
      }
    }
  }


  void calcFluxes() {
    interpForcingInTime();
    super.calcFluxes();
  }


  void updatePrivateParams() {
    sinkingRate[getVarIndex("SD")] = getParamVal("wsink_SD");
    sinkingRate[getVarIndex("LD")] = getParamVal("wsink_LD");
    super.updatePrivateParams();
  }
  

  float attTot(int k, int m) {
    float att = getParamVal("attsw");
    att += getParamVal("attP") * (getVar("SP").current[km2ind(k,m)] + getVar("LP").current[km2ind(k,m)]);
    return att;
  }
  

}






class ForcingSet {

  float[] forcingTime, forcingZ;
  float[][] PARforcing, forcingH, iceForcing;
  float[][][] Tforcing, kappaForcing;
  
  ForcingSet(String forcingFilename) {
    NetcdfFile nc = nc_open(forcingFilename);
    forcingTime = nc_read1D(nc,"time");
    forcingZ = nc_read1D(nc,"z");
    PARforcing = nc_read2D(nc,"par");
    forcingH = nc_read2D(nc,"H");
    Tforcing = nc_read3D(nc,"temp");
    kappaForcing = nc_read3D(nc,"kv");
    iceForcing = nc_read2D(nc,"icecover");
    nc_close(nc);
  }
  
  int length() {return forcingH.length;}
  
  float[] time() {return forcingTime;}
  float[] z() {return forcingZ;}
  float[] one_H_case(int i) {return forcingH[i];}
  float[] one_PAR_case(int i) {return PARforcing[i];}
  float[][] one_T_case(int i) {return Tforcing[i];}
  float[][] one_kappa_case(int i) {return kappaForcing[i];}
  float[] one_ice_case(int i) {return iceForcing[i];}
  
}
